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Abstract 

Motivated by a recent paper by the Potsdam numerical relativity group, we 
have constructed a new numerical code for hydrodynamic simulation of axisym- 
metric systems in full general relativity. In this code, we solve the Einstein field 
equation using Cartesian coordinates with appropriate boundary conditions. On 
the other hand, the hydrodynamic equations are solved in cylindrical coordinates. 
Using this code, we perform simulations to study axisymmetric collapse of rotat- 
ing stars, which thereby become black holes or new compact stars, in full general 
relativity. To investigate the effects of rotation on the criterion for prompt col- 
■ lapse to black holes, we first adopt a polytropic equation of state, P = Kp r , 

where P, p, and K are the pressure, rest mass density, and polytropic constant, 
with r = 2. In this case, the collapse is adiabatic (i.e., no change in entropy), and 
we can focus on the bare effect of rotation. As the initial conditions, we prepare 
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rigidly and differentially rotating stars in equilibrium and then decrease the pres- 
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sure to induce collapse. In this paper, we consider cases in which q = J/Mg < 1, 
£h ■ where J and M g are the angular momentum and the gravitational mass. It is 

found that the criterion of black hole formation is strongly dependent on the an- 
gular momentum parameter q. For q < 0.5, the criterion is not strongly sensitive 
to q; more precisely, if the rest mass is slightly larger than the maximum allowed 
value of spherical stars, a black hole is formed. However, for q <^ 1, it changes 
significantly: For q ~ 0.9, the maximum allowed rest mass becomes ~ 70-80% 
larger than that for spherical stars. These findings depend only weakly on the 
rotational profiles given initially. We then report the results for simulations em- 
ploying a r-l&w equation of state P = (T — l)pe, where e is the specific internal 
energy, to study effects of shock heating. We find that the effects of shock heating- 
are particularly important for preventing prompt collapse to black holes in the 
case of large q [i.e., q = 0(1)]. 



§1. Introduction 

In the 1980s, one of the most important issues in numerical relativity was to per- 
form simulations of rotating stellar collapse with the assumption of axial symmetry. 
It seems that there were two motivations for such studies. cP One was to explore the 
final fate after gravitational collapse of rotating stars. If the cosmic censorship con- 
jecture suggested by PenrosecP is correct, any collapsed star should be surrounded 
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by an event horizon, and consequently settle <down to ja rotating Kexr black hole, 
according to the uniqueness theorem of Israel, 0) Carter a* and others, a' However, if 
this conjecture is not correct, such a star might be reduced to a state with naked 
singularities. Numerical relativity was required to resolve this problem. The other 
motivation was to compute gravitational wave forms in the collapse, because the 
gravitational collapse of a rotating star to a black hole or a neutron star is one of 
promising sources of gravitational waves. 

Simulations of rotating stellar collapse in full general relativity were first per- 
formed by Nakamura. EP Using the (2+l)+l formalism developed by Maeda et al.,0P 
he performed simulations of rotating collapse for massive stars whose masses are 
much larger than the maximum allowed mass for neutron star formation. He used 
cylindrical coordinates (w, z) with a non-uniform grid spacing. At that time, he was 
able to use at most a (42,42) grid resolution for (zu, z) because of restricted compu- 
tational resources. The interesting finding in his simulations was that the rotational 
parameter q = J/Mjj, where J and M g are the angular momentum and gravitational 
mass, can be an important parameter for determining black hole formation: The 
results of his simulations suggest that if q > 1, no black hole is formed, indicating 
that cosmic censorship holds. Unfortunately, he was not able to compute gravita- 
tional waves, because his computational resources were severely restricted and the 
formalism he used was not well suited for computing gravitational wave forms in the 
wave zone. 

To compute gravitational waves emitted during gravitational collapse to black 
holes, Stark and PiranEP performed simulations similar to those of Nakamura, adopt- 
ing spherical polar coordinates with a typical grid size (100, 16) for (r, 9). The distin- 
guishing feature of their work is that they adopted the Bardeen-Piran formalism, W 
which is well suited for computation of gravitational waves in the wave zone. As a 
result of this choice of formalism, they succeeded in computing gravitational wave 
forms, and clarified that the wave forms are characterized by the quasi-normal mode 
of rotating black holes formed after gravitational collapse and that the total radi- 
ated energy of gravitational waves is at most 0.1% of the gravitational mass of the 
system. %> 

Since the completion of their work, no new work in this field has been done 
for about 15 years. Although several questions they originally wished to answer 
have been answered by their simulations, there are still many unsolved issues in 
astrophysics and general relativity that can be investigated using axisymmetric hy- 
drodynamic simulations in full general relativity. For example, the following issues 
should be addressed as open questions: (1) Previous simulations were performed 
using restricted initial conditions and equations of state. tf'EP It is not clear whether 
their results are independent of the initial velocity field, density profile, internal en- 
ergy distribution, and equations of state. A systematic study considering different 
initial settings should be carried out for detailed exploration of the validity of the 
cosmic censorship conjecture. (2) Previous works have indicated that if q > 1, a 
black hole is not formed. However, these works did not carefully investigate how 
the criterion of black hole formation changes due to the effect of rotation for q < 1. 
(3) Stark and Piran computed gravitational wave forms in the case of black hole 
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formation. There have been also many Newtonian works for computing gravita- 
tional wave forms in the formation of neutron stars that are not very massive and 
compact. IHP However, there has been no fully general relativistic work for computing 
gravitational wave forms in neutron star formation. For the case in which the formed 
neutron star is so massive and compact that general relativistic effects play an im- 
portant role even at the formation, fully general relativistic studies on gravitational 
waves could contribute to the understandings of gravitational-wave-astronomy. (4) 
Realistic simulations of the rotating collapse of massive stars, which thereby become 
black holes (or proto-neutron stars) in full general relativity, have not yet been per- 
formed. This should be done to understand mechanisms of black hole formation in 
nature. Actually, a study of the formation of rapidly rotating black holes with sur- 
rounding accretion disks in stellar core collapse is currently one of the hot topics in 
connection with a hypothetical scenario for the central engine of 7-ray bursts. (5) 
The formation mechanism of supermassive black holes of mass 10 6 -10 9 Mq, where 
M & denotes the solar mass, is also an open issue. Recently, Baumgarte and Shapiro 
have proposed an interesting scenario in which rigidly rotating supermassive stars 
at mass shedding limits (i.e., rigidly rotating stars of maximum rotational velocity) 
are possible progenitors of supermassive black holes. According to their scenario, 
after dissipating the internal energy due to thermal radiation, the supermassive stars 
become unstable at the point at which R/M g ~ 450 and q ~ 0.97, where R denotes 
the circumferential radius at the equator. Since q ~ 1, it is not clear whether such a 
progenitor really collapses to a black hole. To confirm that the scenario is promising, 
simulation in full general relativity is necessary. 

All these arguments obviously remind us that there are still many astrophysical 
and general relativistic motivations for performing hydrodynamic simulations of the 
rotating collapse of massive and supermassive stars in full general relativity. 

Before proceeding, we roughly summarize the present situation regarding our 
computational resources. Recently, we have been performing 3D numerical sim- 
ulations in full general relativity with Cartesian coprilinates (x, y, z) for investi- 
gating merging processes of binary neutron stars. tS' ll^P Our 3D simulations are 
performed on the FACOM VPP/300R at the National Astronomical Observatory 
of Japan (NAOJ), in which we can use about 30 GBytes memory. We assume a 
plane symmetry with respect to the equatorial plane [z = 0) and use a grid of size 
(2N + 1, 2N + 1, A + 1) in (x, y, z) where N is at most 160. Let us here estimate a 
desirable grid number N for simulation of the merging of binary neutron stars. In 
the following, we assume to use a uniform grid. 

First, we wish to locate the outer boundary in the wave zone to compute the 
gravitational wave forms accurately. The typical wavelength of gravitational waves 
just before the merging of two neutron stars is about 50M g . Thus, the outer bound- 
aries should be located for > 50M g . Second, such a binary could result in a black 
hole, and to resolve the black hole formation accurately, we need at least 10 grid 
points in the gravitational radius ~ M g . These considerations imply that if we use a 
uniform grid, we need N > 500 for an accurate computation. In the next generation 
machines, which will be available in a year, we will be able to use several hundred 
GBytes. In this case, N = 400(Memory/400GBytes) 1//3 would be feasible. This is 
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close to the value needed for accurate computation. However, with this maximum 
number, it would take ~ 100 CPU hours for one simulation on such machines, in 
which the computational speed is several hundred GFlops. This implies that a sys- 
tematic survey, performing simulations for <^ 100 models of different equations of 
state and compactness of neutron stars, would still be a difficult task to complete on 
a time scale of approximately a year. (Note that under normal circumstances, we 
can use at most ~ 1000 CPU hours per year.) 

In contrast, the memory and computational speed of present computational re- 
sources are large enough for axisymmetric simulations with ~ (500,500). In the 
following, we assume the use of cylindrical coordinates (w,z), with reflection sym- 
metry with respect to z = 0; i.e., the computational domain is assumed to be in the 
region satisfying < w < L x and < z < L z , where L x and L z are constants. In 
this simulation with grid size ~ (500, 500) can be performed on a machine 

with <, 1 GBytes and ~ 1 GFlops using only ~ 10 CPU hours or less (see below). 
This implies that a systematic survey performing simulations for ~ 100 models with 
N ~ 500 can be completed for axisymmetric systems in several months. Using more 
powerful machines, such as FACOM VPP/300R, we could use grids of size about 
(3000,3000). In this case, completion of highly resolved simulations would also be 
feasible. 

One longstanding problem for axisymmetric simulations in full general relativity 
has been to develop methods in which the accuracy and stability for a long-duration 
simulation can be preserved. In axisymmetric simulations, we have in general used 
cylindrical and spherical coordinate systems, which have coordinate singularities at 
the origin and along the symmetric axis w = 0. At such singularities, we have often 
faced several troubles. For example, there appear singular terms in the (3+1) form 



where jij denotes the three metric. Because of the regularity condition along the 
symmetric axis w = 0, these terms are finite in reality, but it is necessary to introduce 
appropriate prescriptions to guarantee them to be finite. Moreover, at the coordinate 
singularities, we have to change the scheme for finite differencing, because there are 
no negative value of w and r. Under such changes, the numerical system can have 
a mismatching between the coordinate singularities and their neighborhoods, often 
resulting in the accumulation of numerical error and, thus, numerical instability. For 
a long-duration simulation, we have often been required to add artificial viscosities 
around the coordinate singularities to stabilize the numerical system.© 

Recently, the Potsdam numerical relativity group has proposed a method to 
make a robust code for solving the Einstein field equation for axisymmetric systems. 
tZP The essence of their idea is that Cartesian coordinates could be used even for 
axisymmetric systems if the Einstein field equation is solved only for y = (or 
x = 0), using the boundary condition at y = ±Ay (or x = ±Ax) provided by 
the axial symmetry. Since the Einstein field equation is written in the Cartesian 
coordinate system, we neither have singular terms nor have to change the finite 
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differencing scheme anywhere, except at the outer boundaries. Thus, it is possible 
to perform a stable and accurate, long-duration simulation without any prescription 
or artificial viscosities. Another merit for us is that it is easy to construct such a 
code, because we have already possessed a 3D code and only need to modify its 
boundary conditions. 

Using this method, we have constructed a new numerical code for axisymmetric 
spacetimes and performed simulations of rotating stellar collapse. As a first step, 
in this paper, we focus on simulations adopting simple initial conditions and simple 
equations of state to investigate the effects of rotation on the criteria for prompt 
collapse to black holes. As initial conditions, we prepare rigidly and differentially 
rotating stars in equilibrium states which are not too compact. We prepare a large 
number of equilibrium states that have different values of q but identical values 
of the rest mass. Then, we reduce the pressure to induce gravitational collapse. 
We pay particular attention to collapse for q < 1 in this paper [i.e., we focus on 
the issue (2) mentioned above] , because black holes are not likely to be formed 
for q > LErlP For most computations, we adopt a polytropic equation of state 
[see Eq. ( |2-2j )]; namely, we do not take into account the effect of shock heating. 
In this treatment, it is easy to investigate the bare effect of rotation for black hole 
formation. For some computations, we adopt a P-law equation of state [see Eq. Q2-3Q 1 
to investigate the effect of shock heating. To roughly mimic a moderately stiff 
equation of state for neutron stars, the adiabatic constant r is set to 2. Giving the 
numerical results for this experiment, we demonstrate that our present strategy and 
current computational resources are sufficient to systematically perform stable and 
well-resolved hydrodynamic simulations in axisymmetric numerical relativity. 

The paper is organized as follows. In §2 we describe the formulation and nu- 
merical methods. In §3 we give the initial conditions which we adopt in this paper. 
In §4 we present the numerical results. Section 5 is devoted to summary and dis- 
cussion. Throughout this paper, we adopt units in which G = 1 = c where G and c 
denote the gravitational constant and speed of light. We use Cartesian coordinates, 
x k = (x, y, z) as the spatial coordinates, with r = \/x 2 + y 2 + z 2 , w = \/ x 2 + y 2 and 
ip = tan _1 (y/x). t denotes the coordinate time. In the following, we use a uniform 
grid with grid spacing Ax = Ay = Az or Ax = Ay = 2Az. Greek indices /i, v, ■ ■ ■ 
denote x, y, z and t, small Latin indices i, j, k, ■ ■ ■ denote x, y and z, and capital Latin 
indices A and B denote x and z. 



§2. Methods 



2.1. Summary of formulation 

For solving the Einstein field equation, we use the same formulation and gauge 



condition as in Refs. 15), 14) and 1£), to which the reader may refer for details and 



basic equations. The fundamental variables for the geometry used in this paper are 
as follows: 



a : lapse function, 
(3 k : shift vector, 
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jij : metric in 3D spatial hyper surf ace, 

7 = e 120 = det( 7 y), 

Ifij = ^ ^Tu'i 
i^ij : extrinsic curvature. 

The variables, 7$,-, 0, Aij = e~^{Kij — jijK k k /3), together with the three auxiliary 
functions Fj = $ dfrjij and the trace of the extrinsic curvature K k k , evolve according 
to a free evolution code. 

As the source of the Einstein equation, we use an ideal fluid. The fundamental 
variables are defined as follows: 

p : rest mass density, 
e : specific internal energy, 
P : pressure, 
h = 1 + e + P/p, 

: four velocity; % = hit}., 



,k 



U 



'•" = -ft; m = v tp , 



w = au 







l + ^muj. (2-1) 



In this paper, we adopt two kinds of equations of state for the dynamical simu- 
lations. One is a polytropic equation of state, 

P = Kp r , (2-2) 

where K is fixed to a constant value throughout the evolution. The other is a /Maw 
equation of state of the form 

P={r-l)pe. (2-3) 

In both cases, we set the adiabatic constant as r = 2 (i.e., the polytropic index n 
is given by n = 1/(7^ — 1) = 1) as a qualitative approximation of moderately stiff 
equations of state for neutron stars. We note that if we prepare a polytropic star as 
the initial conditions and the system evolves in an adiabatic manner without shocks, 
the equation of state is preserved in the form of Eq. (2-2), even using Eq. ( p-3[ ); i.e., 
the value P/p r = K'(x^) for any fluid element remains a constant (= K). 



In the simulation with the equation of state ( |2-2| ), the entropy of each fluid 
element is conserved during the time evolution. For this reason, we refer to the 
collapse with Eq. ( |2-2j ) as "adiabatic collapse" in the following. 

The slicing and spatial gauge conditions we use in r Lhis_paper are basically the 
same as those adopted in our previous series of papers; EJ)'lliP'll3) i.e., we impose an 
"approximate" maximal slice condition (K k k ~ 0) and an "approximate" minimum 
distortion (AMD) gauge condition (D^dt^) — 0, where Di is the covariant deriva- 
tive with respect to ! y%j). For the case in which a rotating star significantly contracts 
and 4> increases by a large factor, we slightly modify the spatial gauge condition in 
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order to improve the spatial resolution around the highly relativistic region, as in 



Refs. [L8D and |19|). In this paper, we modify the gauge according to 



/? x = /3amd-/W^/?amd, (2-4) 
w + e 

P^PiuD-fV^r/AMD, (2-5) 

^ = /3amd[1 -/(*)], (2-6) 

where e is a constant much less than Ax and Az, /?amd ^ s @ k m the AMD gauge 
condition, /3^d = (f3% MD x + (3 y AMD y)/(m + e), and f(t) is chosen as 

1 for 0o > 1.0, 

f(t) = { 200-1 for 0.5 < O < 1.0, (2-7) 
for 0o < 0.5, 

where 0o is at r = 0. Namely, for large 0o > 1, f3 z = and f3 m ~ 0, but f3 v ^ 0. 

We impose axial symmetry, following Ref. |l7|). First, we define the computa- 
tional domain as < x, z < L and — Ay < y < Ay, where L denotes the location 
of the outer boundaries, and reflection symmetry with respect to the z = plane is 
assumed. In this computational domain, we need only 3 points in the y direction, 
and ±Ay. We here determine that the Einstein equation is solved only in the y = 
plane. Then, the boundary conditions at y = ±Ay that are necessary in evaluating 
y-derivatives are supplied from the assumption of axial symmetry as 

Qxx= Q£°J cos 2 cp - Qg) sin 2cp + Qg> sin 2 p, 

Qxy= \Q® sin 2p + Qg) cos 2<p - sin 2p, 

Qyy= Qi°i sin 2 tp + Q®j sin 2p + QjJ cos 2 

Qk= Qi°2 cos 9? - Q{J sin p, Q yz = Q<£] sin tp + cos ip, 
Q x = Q(. 0) cos - Q y 0) sin 99, Q y = sin + Q y 0) cos 99, 

Qzz = Qzz 1 Qz = Qz > Q = 1 

where tp = tp(x) = tan -1 [± Ay / \/x 2 + (Ay)' 2 ], Qij, Qi and Q are an arbitrary tensor, 
vector and scalar in the 3D spatial hypersurface, and Q^' , and denote 
the values of Qij, and Q at (y 'x 2 + (^y) 2 , 0, z), which are interpolated using 
Lagrange's formulac2P with three nearby points along the x direction (i.e., x db Z\x 
and x). At x = L, we use only two points, x — Ax and x, for the extrapolation. We 
note that in using the interpolation and extrapolation, we assume that the geometric 
variables are smooth functions. This assumption is justified as long as black holes 
are not formed. In the case of black hole formation, however, some of the geometric 
variables around the black hole horizons become very steep. In such a situation, the 
interpolation (and extrapolation) becomes less accurate. 

To impose the gauge conditions, as well as to solve constraint equations in 
preparing the initial conditions, we solve scalar and vector elliptic-type equations 



of the form0) 

AQ = S, (2-8) 

AQi = S u (2-9) 

where A denotes the Laplacian in the flat 3D space, and S and Si denote the source 

terms. Using the interpolation mentioned above, d yy Q and d yy Qi are evaluated in 
the finite differencing as 

Q(°) -Q(x,0,z) a „ O z 0) -Q z (x,0,z) 



dyyQ -2 (z ^ )2 , d yy Q z -* ^ 2 

„ n _ QaP\ cosy | - Q x (x,0,z) _ ^ Q y 0) | cosy | - Q y (x,0,z) 

Oyy^ x - z {Ay) 2 ' yy ^ y ~ (Ay) 2 

On the other hand, the finite differencing in the x and z directions, d xx Qi and d zz Qi 
are written in the standard form as 

Qi(x + Ax, 0, z) - 2Qi(x, 0, z) + Qi(x - Ax, 0, z) 



(Ax) 2 

Qi(x, 0,z + Az) - 2Qj(x, 0, z) + Qj(x, 0,z- Az) 
(Az) 2 

Thus, in the finite differencing form for each component of Eq. (2-9), only one com- 
ponent of Qi is included, implying that each component of the vector elliptic-type 
equation is solved independently, as in the case of the scalar elliptic equation. 

For computation of hydrodynamic equations, on the other hand, we do not use 
the above method. Instead, we solve the equations in cylindrical coordinates. These 
equations are written in the y = plane as 

d t p* + ~d w (p*xv x ) + d z (p,v z ) = 0, (2-10) 
x 

dt(p*UA) + -d x (p*xuAV x ) + d z (p*iiAV z ) 
x 

= w!^_ ae 6^ (P + Part) 



ae ^UiUj „ 2ah(w 2 — 1) 

r d A 7 3 

zwn w 

1 



whd A a - ujd A ^ + „,r J d A f j - — V — ~ ' dA4> 



(2-11) 

d t (p*Uy) + ^d x (p*x 2 u y v x ) + d z (p*u y v z ) = 0, (2-12) 



.•><' 1 " " ' 

x A 

d t e* + -d x (e :t xv x ) + d z (e*v z ) = e art , (2-13) 
x 



where = pwe^ and e* = (pe) x l r we°^ . P ar t and e ar t denote artificial viscous 
terms, and are added in the same manner as those in Ref. |l4|) only for the F- 
law equations of state to consider effects of shock heating. Equations (2-10)-(2T3) 
represent the continuity equation, Euler equation for the x and z directions (A = x or 
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z), Euler equation for the y direction (or angular momentum equation), and energy 
(or entropy) equation, respectively. The hvdrodynamic equations are solved using 
van Leer's scheme for the advection termsEH' in this paper. We note that for the 
polytropic equations of state ( |2-2| ), we do not have to solve Eq. (2- 13) B2P nor add the 
artificial viscous terms. 

There are two reasons why we adopt these equations instead of the equations in 
Cartesian coordinates. One is that in the simulation, conservation of the rest mass 
and angular momentum of the system during the evolution is easily guaranteed in 
these forms with help of standard numerical methods. The other is that the hydro- 
dynamic variables often become non-smooth in the presence of shocks in contrast to 
the geometric variables. In such cases, the axisymmetric conditions imposed by the 
interpolation (and extrapolation) are less accurate. 

With these axisymmetric equations, however, there exist coordinate singularities 
along the symmetric axis (x = 0) which might cause a numerical instability for a long- 
duration simulation. In order to avoid such instability, we here add small artificial 
viscous terms 

v v i S p*Au A (2-14) 

to Eq. (2-11) in order to stabilize the system. Here, A denotes the flat Laplacian 
again, and i<Vi s is a small constant such that the typical dissipation time scale due 
to the viscosity, (Ax) 2 /u v i S , is much longer than the dynamical time scale of the 

—1/2 

system, ~ p c ' , where p c denotes the central density at t = 0. 

We note that there are no y-derivatives in Eqs. (2-10) — (2T3), so that we do not 
have to assign grid points at ±Ay for the hydrodynamic variables (i.e., we can save 
memory) . 

To check the accuracy of the numerical results, we monitor the conservation of 
the total rest mass M*, gravitational mass M g , and angular momentum J, which are 
computed in the y = plane as 

M* = An / xdx / dzp*, (2T5) 



o Jo 



,(2-16) 



M g = -2 xdx / dz -2irEe^ + —R <^ K^K^ - {K k k f \ 

Jo Jo 8 8 L J 

J = Att x 2 dx / dzp^iiy, (2T7) 
Jo Jo 

where E = phw 2 — P and R is the Ricci scalar with respect to jij. M* should be 
conserved in all systems. Because of the axial symmetry, J should also be conserved. 
On the other hand, M g should not be conserved because of gravitational radiation. 
However, the total radiated energy of gravitational waves is likely to be quite small 
in the axisymmetric spacetime, so that we can consider M g as an approximately con- 
served quantity. We also monitor the conservation of the specific angular momentum 



spectrum t 



M*(j'o) = 47r/ xdxdzp*, (2-18) 
J j<jo 
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where j is the specific angular momentum computed CIS XUy (= hu v ) and jo denotes 
a particular value for j. We stop simulations before violation of the conservation of 
these quantities is large. 

2.2. Test simulations 



Following Ref. 14), we have performed several test simulations, such as the 
spherical collapse of a pressure-free dust to a black hole, radial oscillation of stable 
spherical stars, and quasi-radial oscillation of stable rotating stars. As in Ref. 14), 



we have confirmed that the collapse of dust can be accurately computed up to the 
formation of the apparent horizon and that the frequency of the fundamental mode 
for radial oscillation of spherical stars is obtained. We have also checked that the 
frequency of the fundamental mode for quasi-radial oscillation of rapidly rotating 



stars agrees with that computed in Refs. 14) and 19). From these results, we are 



convinced that our axisymmetric code is sufficiently robust to perform long-duration 
simulations stably and to provide reliable results. 

§3. Initial conditions for rotating stars 

Initial conditions are prepared in the following manner. First, we adopt rigidly 
as well as differentially rotating stars in (approximately) equilibrium states using 
the polytropic equation of state P = K[p r , with T = 2 and K\ = \. Approximate 
equilibrium states are obtained by choosing a conformally flat spatial metric, i.e., 
assuming 7^ = e^dij. (See, e.g., Ref. |l4|) for the equations to be solved.) This 
approach is computationally convenient, and, as illustrated in Ref. |4|), it provides an 
excellent approximation to exact equilibrium configurations. In particular, we here 
prepare equilibrium states which are not very compact, so that effects of the neglected 
spatial components of the metric are not important. Hence, the approximate initial 
conditions can be almost exact equilibrium states. We also note that this approach 
provides valid initial data for full general relativistic simulations in the sense that 
the solution satisfies the constraint equations exactly. 

Following previous studies, c3 we fix a differentially rotational profile according 

to 

F(f2) = u°u v = A 2 {Q - (2), (3-1) 

where A is an arbitrary constant which describes the length scale over which Q 
changes, and Qq is the angular velocity on the rotational axis. In this paper, A is 
chosen to be w e /3 and w e , where w e denotes the coordinate radius of rotating stars 

"In 

reduces to 

A 2 



at the equator. In the Newtonian limit u° — > 1 and — ► zu 2 f2, the rotational profile 



a = ^o-^-ro • (3-2) 



UJ 



+ A 2 



Thus, for A = va e /3 and w e , O at w = w e (hereafter referred to as £2 e ) is about one 
tenth and half of i?o> respectively (cf., Tables II and III). Note that for A — > 00, the 
rotational profile approaches the rigid rotation. 

We prepare a large number of equilibrium states, changing the angular momen- 
tum parameter of the system, q, but fixing the rest mass M* to 0.05. In Tables I— III, 
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several quantities which characterize the properties of rotating stars are listed for 
A = oo (rigid rotation), A = w e /3, and A = w e . Here, the rotational kinetic energy 
T and the gravitational binding energy W(> 0) are defined as 

T= ~ J P *jVd 3 x, (3-3) 
W= J p*ed 3 x + T + M*-M g . (3-4) 

We note that for unit in which K\ = 1, the maximum allowed rest masses of spher- 
ical stars (hereafter M* smax ) and rigidly rotating stars are about 0.180 and 0.207, 
respectively, with R/M g ~ 5.c3 Hence, the rotating stars prepared here as initial 
conditions are not very compact and massive. We also note that for a rigidly rotating 
star at the mass shedding limit (i.e., a rigidly rotating star of maximum rotational 
velocity), the ratio of the polar coordinate length to the equatorial coordinate length 
is -Rratio — 0.57 for M* = 0.05 with q ~ 1.22. For differentially rotating stars with 
A = w e /3 and w e , even a toroidal star of -Rratio = and q S> 1 can be obtained. 
Since -R ra tio > 0.65 and q < 1, it is found that the rotation of stars here is not very 
rapid. 

We note that T/W for differentially rotating stars is larger than that for rigidly 
rotating stars for an identical value of q. During collapse, T/W is expected to 
roughly increase in proportion to R~ l when J and M g are approximately conserved. 
A star for which T/W is larger than the critical value (~ 0.25) will be unstable 
against the non-axisymmetric dynamical instability. EZP Thus, even if q and M* are 
identical, a collapsing star of differentially rotating initial conditions is likely to be 
more susceptible to the non-axisymmetric instability than a rigidly rotating one. 

To induce gravitational collapse, we initially reduce the pressure by changing K 
from K\ = 1 to a value less than unity without changing the profile of p* and U{. 
(Note that we recompute the constraint equations whenever we modify the initial 
equilibrium configurations.) In the case of the polytropic equations of state with 
units for which G = 1 = c, M* = M*K~ n / 2 is a non-dimensional quantity, and the 
maximum allowed value is 0.180 (= M* S max) for spherical stars and 0.207 for rigidly 
rotating stars, as mentioned above. Thus, if K~ l l 2 is larger than about 0.18/0.05, 
a collapsing star oLsmall q (i.e., q <C 1) is likely to collapse to a black hole in the 
absence of shocks.^ For rapidly rotating stars with q = O(l), the criterion for 
K~ l l 2 is expected to be larger by O(l), due to effects of rotation, because for q > 1, 
no black hole would be formed. H'U 1 In the presence of shocks, K'(= P/p r ) increases 
at shocks, so that the threshold value of K~ 1 / 2 for black hole formation is increased. 

We should mention that we do not consider this numerical experiment to be 
realistic for description of stellar collapses in nature. In realistic stellar collapses, the 
equations of state are much more complicated, effects of cooling and heating due to 
microscopic processes are important, the initial stellar radius is much larger, and the 
initial density configuration and rotational profile may be much more complicated. 
In this paper, we focus on the dynamics of rotating collapse in general relativity, and 
pay particular attention to qualitative effects of rotation on the criterion for prompt 
collapse to black holes. We expect that the present study with its simple equations 
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of states and initial conditions will be helpful for understanding qualitative features 
of realistic rotating stellar collapses. 

§4. Numerical results 

In numerical simulations, we basically adopt a fixed uniform grid in which Ax = 
Ay = Az =constant with grid size (x,z) = (N + 1, N + 1). For collapse of rapidly 
rotating stars, the products after the collapses are often accompanied by very thin 
disks. To improve resolution in the direction perpendicular to the disk, we choose 
Ax = Ay = 2Az =constant using grid size (N + l,2iV + 1) in such cases. To 
check the convergence of our numerical results, we adopt different grid resolutions 
as N = 180, 270, 360 and 500. Typically, we choose N = 360 as the maximum 
grid number. Irrespective of the grid size, the grid covers the same computational 
domain, < x, z < L; namely, the resolution is improved by increasing N. We 
always prepare rotating stars whose surface at the equator {w = w e ) is located at 
2L/3 as initial conditions. 

In numerical simulations, we have used the FACOM VX/4R machine in the data 
processing center of NAOJ. This is a vector-parallel machine composed of four vector 
processors, but we have used only one of four processors in the present numerical 
experiment. The typical memory and CPU time in one simulation with typical grid 
size (361,361) are about 400 MBytes and 25 hours for 15000 time steps. (We note 
that with help of parallelization, the CPU hours could be reduced, although we have 
not yet carried this out for the present numerical code.) 

4.1. Adiabatic collapse for rigidly rotating initial data 

Simulations of adiabatic collapses using polytropic equations of state and rigidly 
rotating initial data sets have been performed for models (B)-(G) described in Table 
I. For each set of initial conditions, the simulations were carried out changing the 
decrease factor of the pressure, K, over a wide range. In Fig. 1, we summarize the 
products after collapse: The horizontal and vertical axes represent q and K^ 1 / 2 , 
respectively. The filled circles and crosses indicate that the products are a new 
star and a black hole, respectively. The filled triangles imply that we were not 
able to judge the nature of the products because a very thin disk is formed in the 
equatorial plane, which cannot be well resolved even with the maximum grid number. 
c3 (As for the process of the disk formation, see, e.g., the first four panels of Figs. 4 
and 5.) The dotted line denotes K~ 1 / 2 = 0.18/0.05. The dashed curve denotes 
the approximate threshold which distinguishes between new star and black hole 
formation. In this numerical experiment, the threshold curve can be approximated 
by a quadratic function of q as K^ 1 / 2 = 1.9g 2 + 3.4 for q < 0.5, and it is a steeper 
function for larger q. 

As we expected (see discussion above), the criterion for K^ 1 ! 2 for prompt col- 
lapse to black holes is ~ 3.6 for small q (i.e., q < 0.5). This implies that for q < 0.5, 
the effect of rotation is not very important for supporting the self-gravity of formed 
rotating stars. However, with increasing q, the critical value of K~ l l 2 increases, and 
for q ~ 0.9, it becomes ~ 6-6.5, which implies that a black hole is not formed after 
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collapse even is ~ 75% larger than M* smax . For model (B), the 

critical value of K~ 1 / 2 can be ~ 7 or even larger; i.e., M* of the resulting star can 
be twice as large as M* S max- For this model, the maximum allowed value of K^ 1 ! 2 
might be much larger, but it is difficult to perform accurate simulations for such a 
large decrease factor, because the thickness of the formed disk is too small to be well 
resolved. 

In Figs. 2 and 3, we display the time evolution of the central density p(r = 0)0 
and a at the origin for model (C) with K~ 1 / 2 = 5.92 and 6.71 and for model (E) 
with K~ 1 / 2 = 3.74 and 3.87. The results shown here are obtained with (361,361) 
(solid curves), (271,271) (dotted curves), (181,181) (dashed curves), and (361,721) 
(dotted-dashed curve). For model (C) with K^ 1 ! 2 = 5.92 and model (E) with 
K~ 1 / 2 = 3.87, the collapse is halted during the evolution, resulting in an oscillating, 
compact star. On the other hand, for model (C) with K^ 1 ! 2 = 6.71 and for model 
(E) with K~ x l 2 = 3.87, the stars collapse to black holes in the dynamical time scale. 
As Figs. 2 and 3 show, these results depend only weakly on the resolution, implying 
that convergence is achieved fairly well. 

In Figs. 4-7, we show snapshots of the contour lines and velocity fields in the 
x-z planes at selected time steps for model (C) with K~ 1 / 2 = 5.92 (Fig. 4) and 6.71 
(Fig. 5), and for model (E) with K' 1 / 2 = 3.74 (Fig. 6) and 3.87 (Fig. 7). All the 
results presented here were obtained with (361,361) grid resolution. For model (C) 
with larger initial decrease of the pressure, the rotating star collapses to a highly 
flattened object in the early stage. When the density increases sufficiently, the central 
region stops contracting, to form a core. On the other hand, the outer region, which 
forms a thin disk, is still collapsing. If the decrease factor of the pressure is not very 
large, the core subsequently expands, sweeping the surrounding disk, and finally an 
oscillating, oblate star is formed. On the other hand, for a sufficiently large decrease 
factor (i.e., K~ x l 2 ^ 6.5), the core cannot expand, due to the insufficient pressure, 
resulting in a black hole after a fraction of matter accretion. In this case, the disk of 
a small fraction of mass appears to be formed around the black hole at r ~ 3-WM g 
(see also Fig. 9 (a) and discussion below). 

We note that even for model (C), the product is not a flattened object, but 
spheroidal one, if K' 1 / 2 <, 4. In such cases, the core radius is larger, and the 
evolution process is similar to that described by Fig. 6. The formation of a highly 
flattened object is a feature in the collapse for large q [i.e., q = 0(1)] and for near 
critical values of K~ x l 2 (i.e., K~ x l 2 ^> 6), in other words, for the highly supramassive 
case in the sense M* S> -M* smax . 

For model (E), in which q is not very large, on the other hand, the collapse 
proceeds in a nearly spherical manner with tiny disks. When the density sufficiently 
increases, a nearly spherical core is formed. Since the decrease factor is smaller than 
that for model (C), the radius of the core is larger. For K~ 1 / 2 = 3.74, the core stops 
contracting, mainly due to the pressure force, subsequently forming an oscillating, 
nearly spherical star. For a larger decrease factor, the collapse cannot be halted, and 
a black hole is finally formed, swallowing almost all the fluid elements. 

The noteworthing difference between the evolution for models (C) and (E) for 
near critical values of K is with regard to the rest mass of the core in the early stage 
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of its formation. For (C), the rest mass is not very large, because the fraction of the 
rest mass in the surrounding disk is fairly large (see the third panel of Fig. 4). This 
implies that the core does not have to support all of the rest mass of the system, and 
consequently the internal energy to support the self-gravity can be small. On the 
other hand, for (E), almost all the fluid elements accrete to the core simultaneously 
(see the third panel of Fig. 6), and the fraction of the rest mass which the core has to 
support is larger. This appears to be part of reason why the critical value of K~ l l 2 
for model (C) can be larger than that for model (E). 

We should note that the product shown in Fig. 4 can have a large T/W, because 
according to a simple scaling law assuming the conservation of the angular momen- 
tum and the mass, we can expect it to roughly increase as ~ (T/W)i(Ri/R), where 
(T/W)i and Ri are T/W and R at t = 0. Since T/W ~ 0.06 at t = (see Table I), 
T/W in a formed oscillating star would be ~ 0.2-0.3, and non-axisymmetric insta- 
bilities should be taken into account. A formed star is evidently secularly unstable 
with respect to gravitational wave emission to form a bar.EiP However, the secular 

— 1/2 

time scale is much longer than the dynamical time scale, ~ p c , so that this would 

— 1/2 

not be relevant even in a non-axisymmetric simulation performed only for ~ 3p c 



A star also might be dynamically unstable to form a bar or spiral arms.EZP Using the 
results found in a previous paper, c3 we can roughly examine the dynamical stability. 
From Fig. 2, the normalized central density p(r = 0) = Kp(r = 0) is found to be 
between 0.05 and 0.15, while the normalized rest mass M* is about 0.3. In Ref. |27| ), 
we have found that if p <^ 0.1, a star with M* 0.3 and T = 2 is dynamically 



unstable (cf., Fig. 2 of Ref. 27)). Thus, such a star is located near the threshold 



of the dynamical stability. Therefore, the process of the collapse might change if 
non-axisymmetric simulations were carried out. 

In Fig. 8, we display the angular velocity Q as a function of w in the equatorial 
plane at selected times for model (C) with K~ l l 2 = 5.92. It is found that the 
star is differentially rotating, although it was rigidly rotating at t = 0. Since the 
formed star is not stationary but in an oscillating state, 17 varies with time and 
becomes a steeper function when the central density is large. As shown in a previous 
paper, E§ the effect of differential rotation can significantly increase the maximum 
allowed rest mass much beyond M* S max- This seems to be the reason why this highly 
supramassive object of rest mass much larger than M* smax did not collapse to a black 
hole in this simulation. 

We, however, note that the differential rotation will be changed to rigid rotation 
due to a magnetic field or viscosity, or that the angular momentum will be dissipated 
by the magnetic dipole radiation for a realistic neutron star formed after the collapse. 
c3 Thus, even if they are formed, any real supramassive stars of differential rotation 
will eventually collapse to black holes due to dissipation processes in the secular time 
scale. 

In Figs. 9 (a) and (b), we plot the fraction of the rest mass inside a coordinate 
radius r, M*(r)/M*, as a function of time (a) for model (C) with K~ x l 2 = 6.71 and 
(b) for model (E) with K~ 1 / 2 = 3.87. In both cases, black holes are formed on the 
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dynamical time scale. Here, M*(r) is denned as 

M*(r) = 4vr f x'dx'dz'p*. (4-1) 

Jr'<r 

The dotted, solid and dashed curves denote the results for (271,271), (361,361) and 
(361, 721) grid resolutions, respectively. As expected from Fig. 7, almost all the 
fluid elements are swallowed inside the apparent horizon for model (E) . On the other 
hand, a small fraction of mass could form a disk around a black hole for model 
(C). Since the simulation crashed soon after the formation of the apparent horizon 
at t ~ 1.15p~ 1//2 , it is not possible to draw a definite conclusion from Fig. 9(a). 
However, Fig. 9(a) appears to suggest that r j^-j5% of the total rest mass might form 
a disk around a black hole for r ~ 3-10M 9 .B The reason for these results can be 
explained from the specific angular momentum distribution at t = as follows. 

In Fig. 10, we show M* (j) /M* as a function of j/M g for models (C) and (E) (note 
that M g here denotes the initial value). For (C), we show it at selected time steps to 
demonstrate that the specific angular momentum spectrum M*(j) is preserved during 
the evolution fairly accurately. In both models, the specific angular momentum j is a 
monotonically increasing function of w, and at the equator, it is ~ 3.6M g for model 
(C) and ~ 1.8M g for model (E). Taking into account that J is a conserved quantity, 
M g is almost conserved, and the disk mass occupies only a small fraction of the total 
mass, we can expect that q of the formed black hole should be approximately equal 
to the initial value. Thus, the specific angular momentum at the innermost stable 
circular orbit in the equatorial plane around the formed Kerr black holes, which is 
the minimum allowed value for a test particle orbiting a Kerr black hole, is ~ 2.1M g 
for model (C) and ~ 2.9M g for (E).E3) It is immediately found that the specific 
angular momentum of no fluid elements is large enough to form a disk around the 
black hole for model (E). For model (C), a fraction of fluid elements of j between 
~ 2.\M g and 3.6M 9 can form a disk, but Fig. 10 shows that the mass fraction is 
only ~ 6%, which is in approximate agreement with the fraction of the disk mass 
found in the numerical simulation. 

4.2. Adiabatic collapse for differentially rotating initial data 

Simulations of adiabatic collapses using polytropic equations of state and differ- 
entially rotating initial data sets have been performed for models (H)-(K) described 
in Table II and for models (L)-(O) in Table III. For each set of initial conditions, 
the simulations were carried out changing the decrease factor of the pressure for a 
wide range. In Figs. 11 and 12, we summarize the products formed after collapse for 
A = zz7 e /3 and A = w e in the same manner as for Fig. 1. Note that the dashed curves 
denote the approximate threshold for black hole formation (i.e., below the dashed 
curves, a stable star is formed) for rigidly rotating initial data sets. It is found that 
the threshold for the differentially rotating initial conditions almost coincides with 
that for the rigidly rotating ones. Therefore, for the moderate range of rotational 
profiles we have investigated (tu e /3 < A < oo), the criterion of black hole formation 
depends only weakly on the initial rotational profile. 

On the other hand, the products after the collapse depend strongly on the initial 
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rotational profile, in particular, for q ~ 0.9: For rigidly rotating initial data, the star 
remains a spheroid in the collapse even if q ~ 0.9. On the other hand, it becomes 
a toroid for differentially rotating initial data with A = tu e /3. To illustrate this 
fact, we display snapshots of the contour lines and velocity fields in the x-z planes 
at selected time steps for model (H) with K' 1 ! 2 = 6.32 in Fig. 13. In this case, the 
star does not form a black hole, but a toroid. As in the case of rigidly rotating initial 
data, fluid elements which are initially located far from w = collapse toward the 
central region in the early phase. On the other hand, those near the symmetric axis 
do not simply collapse toward the center, but slightly expand in the w direction. 
This seems to be due to the effect of strong centrifugal force near the z axis (note 
that Q at w = is initially about 4 times larger than that for the rigidly rotating 
cases of identical q). As a result of these types of motion, the collapsed object forms 
a ring-shape core in the early stage (cf., the third panel). The ring subsequently 
collapses toward the center, to become a very compact toroid (cf., the fifth and sixth 
panels). However, the centrifugal force and/or pressure hang up the collapse, and 
subsequently the ring expands with the formation of a jet along the symmetric axis 

— 1/2 

(cf., the seventh panel). Once it reaches a maximum expansion at t ~ 1.5p m ax , the 
ring begins to contract again to form a compact toroid (cf., the eighth panel), and 
this oscillating motion is repeated. 

Here, we note that the maximum normalized density of the ring is very small 

— 1/2 

Pmax ~ 0.03 at t ~ 1.5p m ax , while M* is very large, ~ 3.1. This implies that the 
ring would be dynamically unstable with respect to formation of a bar and/or spiral 
arms, or with respect to fission in a non-axisymmetric simulation. EZP This result is 
different from that in the case of rigidly rotating initial data, in which a collapsed 
object seems to be marginally stable, and it suggests that a collapsed object of 
differentially rotating initial data could be more susceptible to a non-axisymmetric, 
dynamical instability than that of rigidly rotating initial data, even if q is identical 
initially. 

For model (H) with K^ 1 / 2 > 7, black holes are formed. For these cases, the 

— 1/2 

evolution processes until t ~ p m ax is essentially the same as that shown in Fig. 13, 
but the star subsequently collapses, without bounce. In this black hole formation 
case, almost all the fluid elements are swallowed inside the black holes quickly. This 
can be understood from the sixth panel of Fig. 13. This result is in contrast to that 
in the black hole formation for model (C). The reason is as follows. In the case of a 
steep differentially rotational profile at t = 0, the specific angular momentum near 
the equator is much smaller than that for rigidly rotating cases. For model (H), the 
specific angular momentum at the equator is 0.94M 9 , which is smaller than the value 
at the innermost stable circular orbit of a Kerr black hole with q = 0.9 and mass 
M g . Consequently, it is unlikely for a disk to be formed. We argue this point in §5 
in detail. 

In Fig. 14, we show snapshots of the contour lines and velocity fields in the 
x-z plane at selected time steps for model (J) with K~ l l 2 = 3.74, in which a black 
hole is not formed. Since q is not very large (~ 0.48) in this case, a ring is not 
formed throughout the entire evolution. However, the effect of the rotation appears 
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to be significant. Compared with the core formed in the simulation of model (E), 
in which q is nearly equal to that for model (J), the shape of the core is much 
more deformed, due to the effect of the differential rotation initially given (compare 
the second panel with the second panel in Fig. 6). Also, a violent non-spherical 
oscillation is found after formation of the new star. These results imply that even if 
q and M* are identical initially, the shape and motion of the products after collapse 
depend strongly on the initial rotational profiles. 

4.3. Non-adiabatic case 

Simulations of non-adiabatic collapses using a /Maw equation of state have been 
performed for models (C)-(E) and (H)-(J) described in Tables I and II. In Figs. 15 
and 16, we summarize the products after collapses for (C)-(E) and (H)-(J) in the 
same manner as for Fig. 1. The dashed curve denotes the approximate threshold 
for black hole formation in the adiabatic collapses. Obviously, the shock heating 
raises the criterion for K~ l l 2 for prompt formation of black holes. It is found that 
the effect is in particular significant for larger q. We note that for model (C), we 
cannot judge the nature of the product for K~ 1 / 2 > 10, because a very thin disk is 
formed in the equatorial plane and cannot be well resolved even with (501, 501) grid 
resolution. For model (H), formation of black holes is not found for K" 1 ! 2 < 10. 
For larger K~ 1 / 2 (i.e., K~ x l 2 S> 10), black holes would be formed for both models, 
but to find this formation, simulations with finer resolution are necessary to resolve 
the thin disk and/or thin ring. 

Even in the presence of the shock heating, the formation process of black holes 
is not modified drastically. In Fig. 17, we show snapshots of the contour lines and 
velocity fields in the x-z planes at selected time steps for model (E) with K~ l l 2 = 
5.00 as an example. Although shocks are generated at the formation of the central 
core (cf., the second panel), they stall due to the subsequent accretion of matter, 
and hence the collapse is never halted, resulting in a black hole on the dynamical 
time scale. 

On the other hand, the process in the case of no black hole formation is con- 
siderably affected by shocks, in particular, for q ~ 0.9. In Figs. 18-21, we show 
snapshots of the contour lines and velocity fields in the x-z planes at selected time 
steps for model (C) with K~ 1 / 2 = 8.37, for model (E) with K~ 1 / 2 = 4.47, for model 
(H) with K~ 1 / 2 = 9.49, and for model (J) with K~ 1 / 2 = 4.00. 

In the case of small q, spheroidal shocks are formed around the surface of the 
central core (cf., the first panel in Figs. 19 and 21). This implies that the shock 
heating is effective only around the spheroidal shell. Indeed, K'{x^) around the 
center of the core remains nearly constant throughout the evolution. Therefore, the 
shock heating is not very effective for increase of the internal energy in the high 
density region. Instead, the energy generated by shocks at the shell is effectively 
used to provide the energy to the surrounding matter, which subsequently expands 
outward. 

On the other hand, the formation mechanism of shocks for q ~ 0.9 and K~ l l 2 M* ^> 
M tsmax is completely different from that for smaller q (see Figs. 18 and 20): In the 
very early phase, the shocks around the small spheroidal core (for model (C)) or 



18 



disk-plus-ring (for model (H)) are formed. Here, the shape of the shocks depends 
strongly on the rotational profile initially given. The reason why the mass fraction 
around the shock forming, high density region is small is that a large fraction of 
mass, which is initially located at large vo, has not yet accreted, due to the centrifu- 
gal force, and the collapse proceeds mainly in the z direction. Since the mass of the 
high density region is small, prompt formation of a black hole does not occur in spite 
of the small internal energy. 

The shocks formed in the early phase subsequently expand in the z direction for 
model (C) and toward the z axis for model (H), because a large amount of matter still 
accretes near the equatorial plane and prevents the shocks from expanding outward. 
The mass fraction near the symmetric axis is not very large, so that this shock 
expansion is not efficiently used for sweeping the accreting matter around the cores. 
The evolution processes after a substantial fraction of matter residing in the outer 
part has accreted to the small cores are different for the two cases. For model (C), 
the shocks formed in the interface between the small core and the surrounding disk 
initiate expansion that sweeps the matter (cf., the fourth and fifth panels in Fig. 18). 
However, the expansion is not very powerful because a part of the internal energy has 
already been used in the first expansion in the z direction and these second shocks 
are not very strong. For model (H), the evolution process is more complicated. 
During the accretion of the matter from the outer part to the core, the toroidal 
core starts contracting toward the z axis without significant shock formation (cf., 
the third panel in Fig. 20); i.e., in this case, shocks which sweep accreting matter 
are not formed. Then, the inner surface of the toroid collides at the z axis (cf., 
the fourth panel), and as a result, a shock is formed along the z axis. This shock 
subsequently expands mainly in the w direction (cf. , the fifth panel) , contributing to 
an increase of the internal energy for a large fraction of mass elements in the toroid. 
The toroid oscillates nearly periodically, so that this process is repeated (cf., the 
sixth-ninth panels). For both models (C) and (H), several different types of shocks 
are formed in different stages incoherently. Consequently, the motion of the shocks 
is less coherent and less efficient in sweeping the accreting matter than for the small 
q case. This result appears to represent an extreme example in which rapid rotation 
prevents the formation of coherent shocks. 

The above described mechanism of shock formation also explains the reason 
that the threshold value of K~ x l 2 is much larger than that in the absence of shocks 
only for the large q case. For small q, shocks are generated near the surface of a 
central dense core, so that the entropy can be increased only in the surface region 
of low density. Consequently, the shock heating is not very effective for increasing 
the internal energy of the central core. On the other hand, for large q with large 
pressure decrease (i.e., K^^M* 3> M* S max), shocks can be generated near high 
density regions, i.e., around the small cores and the disks. Indeed, K' is increased 
by a large factor for most fluid elements. This is due to the highly non-spherical 
motion induced by the large centrifugal force, as mentioned above. As a result, a 
large amount of the matter is affected by the shock heating, resulting in a large 
increase of the total internal energy. Thus, although rapid rotation prevents the 
formation of coherent shocks, it could contribute to a shock heating that is efficient 
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for preventing prompt collapse to black holes. 

§5. Summary 

We constructed a new numerical code for axisymmetric hydrodynamic simula- 
tions in full general relativity. In this code, the Einstein field equation is solved 
using Cartesian coordinates with appropriate boundary conditions, while the hy- 
drodynamic equations are solved in cylindrical coordinates. It was shown from nu- 
merical experiments that this method is sufficiently robust to perform stable and 
well-resolved computations for rotating stellar collapse. 

Using this code, we have investigated the criterion for black hole formation 
after rotating stellar collapses. It is found that the effect of rotation could play an 
important role in preventing prompt formation of black holes. For the case in which 
the angular momentum parameter q is smaller than ~ 0.5, a black hole is formed if 
the total rest mass is slightly larger than the maximum allowed rest mass of spherical 
stars (M* S max)) implying that the effect of rotation is not very important. On the 
other hand, the threshold of the total rest mass can be significantly increased by the 
effect of rotation for q <^ 1: e.g., for q ~ 0.9, the threshold value is about 70-80% 
larger than M* smax - The self- gravity of such large mass products is supported by 
rapid, differential rotation. These results are found to depend only weakly on the 
initial rotational profile within the range that we have investigated. 

We have found that the formation mechanism of shocks depends strongly on the 
rotational parameter q. For collapse of slowly rotating stars of q < 0.5, shocks are 
formed at the spheroidal shell around a high density core. As a result, shock heating 
is effective only around the shell and is not very helpful for increasing the internal 
energy of the central core. On the other hand, shocks can be formed in a highly 
non-spherical manner and near the high density region in the collapse of rapidly 
rotating, supramassive stars for which q = 0(1) and K~ l / 2 M* S> M* S max- Since the 
centrifugal force is large in this case, a large amount of matter at large w cannot 
accrete immediately, and the high density regions formed in the early stage could 
have highly non-spherical shapes. As a result, shocks can be formed in a highly non- 
spherical manner around the high density region. (The shapes of the shocks depend 
strongly on the initial rotational profile.) Because of the matter accretion around 
the equatorial plane, shocks cannot expand coherently, and hence they cannot sweep 
the accreting matter effectively. However, shocks can contribute to the heating for 
a large number of the mass elements, because they are formed near the high density 
regions. Consequently, they prevent prompt collapse to black holes even for highly 
supramassive stars with K~ l l 2 M* 3> M* S max- 

The numerical results reported in this paper indicate that the disk mass around 
formed black holes is not very large even for q ~ 0.9. In a different setting, with 
different equations of state, different initial radii, and different initial rotational 
profiles, results could be modified. It is an interesting issue to clarify what kind of 
initial conditions are preferable for obtaining a system of a black hoM surrounded 
by massive disks, because the formation of disks around a black hole £3) could play 
an important role in nature. Here, we carry out a simple analysis to address that 
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equations of state and initial velocity profiles of a star before collapse play a crucial 
role in determining the disk mass. 

In the following, we assume that (1) rotating stars at t = are not very compact, 
(2) black holes are formed only for q < 1, and (3) the fraction of the disk mass is 
much smaller than the total mass and that the value of q of a formed black hole is 
nearly equal to the initial value. Due to the assumption (1), the angular momentum 
can be approximately computed from initial data as 

J — J dVp*w 2 Q, (5-1) 

where we use h ~ 1 and j ~ w 2 Q. For simplicity, here, we consider Q to be of the 
same form as in Eq. ([3-2]). For the case that A <C w e , J ~ M*A 2 i?o — M*zu 2 [2 e , 
where Q e is the angular velocity at vo = w e . Then, q can be approximately written 
as m 2 f2 e /M g , where we use M g ~ M*, due to the assumption (1). For the case 
A ^> w e , on the other hand, the star is in a state of nearly rigid rotation (l?o = ^e), 
and J = M*(vu 2 )f2 e , where 

(m 2 ) =^-J p*m 2 dV. (5-2) 

Consequently, q ~ (m 2 )Q e /M g . Then, the specific angular momentum j can be 
approximately written as 

.„2n f M g q for A « w e , 

Thus, j(w)/M g = q for A <C ro e , and j(w)/M g < qw 2 /{w 2 ) for A S> tu e . To form a 
disk around a black hole, j/M p has to be larger than £isco> where ^isco denotes j/M g 
at the innermost circular orbit around a black hole which depends on q. (^isco = 1 
for q = 1, ~ 2.1 for g = 0.9, and ~ 2\/3 for g <C 1.) This implies that for A <C ro e , a 
disk cannot be formed, irrespective of the value of q for q < 1. For a star which was 
initially in rigid rotation, a disk can be formed if qw 2 /(zu 2 ) = 0(1). For equilibrium 
stars with stiff equations of state, the density profile is not very centrally peaked, so 
that w 2 /(w 2 ) is of O(l). This implies that q has to be of 0(1) for disk formation. 
We presented such an example in §4.1. On the other hand, q may be smaller for 
equilibrium stars of soft equations of state for which the density profile is centrally 
peaked, and w 2 /(w 2 ) can be much larger than unity. We note that w 2 /(w 2 ) depends 
weakly on w e for non-relativistic stars of an identical equation of state. Thus, this 
conclusion is almost independent of the initial stellar radius. From these arguments, 
we can recognize that to obtain disks around a black hole, irrespective of the initial 
stellar radius, (i) stars before collapse should not have a steep differential rotational 
profile, but rather a nearly rigidly rotational profile, (ii) the value of q of stars before 
collapse should be of 0(1) if the equation of state is stiff, or stars before collapse 
should have a soft equation of state if q -C 1. In realistic situations, a progenitor of 
a supernova just before collapse has a soft equation of state with r ~ 1.3. Thus, a 
disk would be formed as long as the initial rotational profile is not too steep. 
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Although these conclusions would be qualitative correct, nothing quantitative 
can be clarified at this time. In particular, we assumed that the disk mass is much 
smaller than the black hole mass in this analysis, but this might be incorrect for the 
case of soft equations of state. This is because the star is centrally condensed in this 
case, and it may be the case that only the central region collapses to form a black 
hole of a small mass fraction. Thus, to compute the disk mass in rotating stellar 
collapse correctly, numerical simulations are obviously necessary. 

In this paper, we have studied rotating stellar collapse using simple equations of 
state and simple initial conditions which are not very realistic. We believe that the 
significant effects of rotation found here are qualitatively correct, but they might not 
quantitatively realistic. To obtain a result applicable to realistic systems, we have 
to perform simulations using more realistic equations of state and initial conditions. 
We plan to perform such simulations in the future. 
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Table I. Central density p c , gravitational mass M g , q = J I M 2 a , T/W, R/M g (where R is the 
circumferential radius at the equator), angular velocity Q, and ratio of polar coordinate radius 
to equatorial radius {R T3 x\o) for rigidly rotating stars. M* is fixed to 0.05, and the polytropic 
constant K is fixed to f . 



Pc 


Mg 


q 


T/W 


R/Mg 


Q 


^ratio 


Model 


0.0f8750 


0.04909 


0.978 


0.0715 


30.04 


0.0892 


0.717 


(B) 


0.0f9316 


0.04907 


0.902 


0.0628 


29.07 


0.0861 


0.750 


(C) 


0.020757 


0.04903 


0.702 


0.0411 


27.02 


0.0748 


0.833 


(D) 


0.022f79 


0.04899 


0.473 


0.0200 


25.40 


0.0557 


0.917 


(E) 


0.022868 


0.04898 


0.327 


0.00986 


24.70 


0.0403 


0.958 


(F) 


0.023408 


0.04897 


0.143 


0.00194 


24.19 


0.0183 


0.992 


(G) 



Table II. The same as Table I, but for differentially rotating stars with A — cc e /3. The rest mass 
is also fixed to 0.05. p max , instead of p c , is given in the first column. In the sixth column, 
and fl e are given. 



pmax 


Mg 


<7 


T/W 


R/Mg 


{o ,n e ) 


-Rratio 


Model 


0.015286 


0.04910 


0.900 


0.0767 


27.51 


(0.3668,0.0334) 


0.650 


(H) 


0.017432 


0.04905 


0.715 


0.0528 


26.31 


(0.3611,0.0329) 


0.750 


(I) 


0.020327 


0.04901 


0.487 


0.0268 


25.14 


(0.2470,0.0223) 


0.867 


(J) 


0.022776 


0.04897 


0.230 


0.00640 


24.32 


(0.1265,0.0114) 


0.967 


(K) 



Table III. The same as Table II, but for differentially rotating stars with A = cc e . The rest mass 
is also fixed to 0.05. 



pmax 


Mg 


q 


T/W 


R/Mg 


(n ,n e ) 


-Rratio 


Model 


0.018646 


0.04907 


0.893 


0.0635 


27.76 


(0.1239,0.0580) 


0.767 


(L) 


0.020153 


0.04904 


0.717 


0.0438 


26.50 


(0.1096,0.0512) 


0.833 


(M) 


0.021913 


0.04900 


0.478 


0.0210 


25.18 


(0.0813,0.0379) 


0.917 


(N) 


0.023066 


0.04897 


0.252 


0.00610 


24.38 


(0.0458,0.0213) 


0.975 


(0) 
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Fig. 1. Summary of the products in adiabatic collapse with rigidly rotating initial data. The 
horizontal and vertical axes denote q = J/M 2 and K" 1 ^ 2 . The filled circles and crosses indicate 
that the products are a new compact star and a black hole, respectively. The filled triangles 
imply that we could not judge the nature of the products. The dotted line and dashed curves 
denote K~ 1 ^ 2 = 18/5 and the approximate threshold for black hole formation, respectively. 
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Fig. 2. Time evolution of p and a at r = for model (C) with K~ 1/2 = 5.92 and 6.71. The 
solid, dotted, dashed, and dotted-dashed curves represent the results for (361,361), (271,271), 
(181, 181) and (361, 721) grid resolutions, respectively. 




Fig. 3. The same as Fig. 2, but for model (E) with A" 1/2 = 3.74 and 3.87. 
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Fig. 4. Snapshots of density contours for p» and the velocity flow for v* in the (x,z) plane for 
model (C) with K^ 1 ^ 2 = 5.92. The contour lines are drawn for p»/p* max = 10~°' 4j for j = 
0, 1,2, • • • , 10, where p„ max is 1.15, 422, 2054, 2705, 360, 264, 876, 315 and 244 times larger 
than p» max at t = 0. The length of the velocity vector (the arrows) is normalized to 0.2c. The 
time and length appear in units of p^ 1 ^ 2 and GM g /c 2 at t = 0, respectively. 
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5. The same as Fig. 4, but for model (C) with K ^ 2 = 6.71. The contour lines are drawn 
for p,/p, max = 1CT ' 4J for j = 0, 1,2, • • • , 10, where p* max is 1.15, 598, 7830, 6.37 x 10 4 , and 
3.89 x 10° times larger than p* max at t = 0. The last panel is the magnification of the 5th 
panel. The thick solid curve in the last two panels indicates the apparent horizon. 
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Fig. 6. The same as Fig. 4, but for model (E) with K 1 ' 2 = 3.74. The contour lines are drawn 
for max = l(r ' 4j for j = 0, 1, 2, • • • , 10, where p* max is 1.12, 724, 952, 475, 126 and 861 

times larger than p„ max at t = 0. 
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7. The same as Fig. 4, but for model (E) with K = 3.87. The contour lines are drawn 
for p,/p, max = lCT ' 4j for j = 0, 1, 2, • ■ • , 10, where p, max is 1.23, 1120 , 2670, 1.78 x 10 4 , and 
1.82 x 10° larger than the initial value. The last panel is the magnification of the 5th panel. 
The thick solid curve in the last two panels indicates the apparent horizon. 
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Fig. 8. Q in units of M g 1 at selected time steps for model (C) with K 1 '" 2 = 5.92. Note that the 
rotating stars are oscillating with time, and at t — 1.74pc 1 '" 2 the radius is smaller than at other 
times. 
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Fig. 9. Evolution of the fraction of the rest mass inside a coordinate radius r, M*(r)/M* (a) for 
model (C) with if" 1/2 = 6.71 and (b) for model (E) with K~ 1/2 = 3.87. The dotted, solid 
and dashed curves represent the results for (271, 271), (361, 361) and (361, 721) grid resolutions, 
respectively. 
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Fig. 10. M*(j)/M, as a function of j/M g for models (C) and (E). The solid, dotted and dashed 
curves for model (C) denote the results at tpl^ 2 = 0, 0.901, and 1.10 for K ~ 1 / 2 = 6.71. 
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;. 11. The same as Fig. 1, but for differentially rotating initial data with A = vj e /3. The dashed 
curve represents the approximate threshold for black hole formation for rigidly rotating initial 
data. 
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Fig. 12. The same as Fig. 11, but for differentially rotating initial data with A = tu e - 
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Fig. 13. The same as Fig. 4, but for model (H) with K 1 ' 2 = 6.32. The contour lines are drawn 
for p,/p, max = l(T ' 4j for j = 0, 1, 2, • • • , 10, where p, max is 1.10, 9.22, 77.7, 99.3, 1781, 4131, 
75.1, 85.3, and 3619 times larger than p» max at t = 0. The time appears in units of p max at 
t = 0. 
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Fig. 14. The same as Fig. 13, but for model (J) with K 1//2 = 3.74. The contour lines are drawn 
for p„/p, max = 10- 4j ' for j = 0, 1,2, • • • , 10, where p, max is 1.06, 61.5, 436, 1554, 610, 63.7, 
41.8, 241, and 36.5 times larger than p* max at t — 0. Note that p c — pmax at t = in this case. 
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Fig. 15. The same as Fig. 1, but for non-adiabatic collapse with rigidly rotating initial data. The 
dashed curve represents the approximate threshold for black hole formation in adiabatic collapse. 
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Fig. 16. The same as Fig. 1, but for non-adiabatic collapse for differentially rotating initial data 
with A — G7e/3. The dashed curve represents the approximate threshold for black hole formation 
in adiabatic collapse. 
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17. The same as Fig. 4, but for non-adiabatic collapse for model (E) with K' 1 ' 2 = 5.00. The 
contour lines are drawn for p* /p* max = 10 -0 ' 4 -' for j = 0, 1, 2, • • • , 10, where p* max is 643, 3814, 
2.65 x 10 4 , and 2.26 x 10 5 times larger than p* max at t = 0. The last panel is the magnification 
of the 4th panel. The thick solid curve in the last two panels indicates the apparent horizon. 




18. The same as Fig. 4, but for non-adiabatic collapse for model (C) with K 1 '" 2 = 8.37. The 
contour lines are drawn for p»/p* max = lCP ' 4 -' for j = 0, 1, 2, • • ■ , 10, where p, max is 1178, 
5015, 1186, 403, 317, 619, 814, and 436 times larger than p, max at t = 0. 



39 



12 ETTFf 

1 iririr 

10 
8 

6 h 



2 pm 



2 



t = 0.757p o " 1 / 2 



»?)? i a ? g i ?U Ada, i 'n 
4 6 8 10 U 




4 6 8 10 12 

X / M 



12 
10 



4 m 



t = 0.905p o " 



t=1.05p; 1 / 2 



12 



a a<- i a ^ ? i <-ij?d 'n 

4 6 8 10 12 



jyr 

10 f- 

B 



4 

2 = 







X / M g 
t=1.46p;'/ 2 



2 4 6 8 10 12 

X / M 




2 4 6 8 10 12 

X / M g 

t=1.76p;>/ 2 




12 
10 



4 r~ 

2 Ef 



t = 2.14p; 1 / 2 







2 



■ iM i ill Wi ,Vfa-,44. 
4 6 8 10 12 



12 
10 
8 

6 f 
4 

2 





t = 2.58p; 1 / 2 



X / M„ 



2 4 6 8 10 12 

X / u 



12 
10 







>H M 1 1 1 1 4 


I i i i>i 

L' ' ' ! 




-0 2c: 






•J ; lijTTv 




,'A; tiii 


,(, ,- us:, , t ,„:f 



2 4 6 8 10 12 

X / M. 



Fig. 19. The same as Fig. 4, but for non-adiabatic collapse for model (E) with K ^ 2 = 4.47. The 
contour lines are drawn for max = lCP 0,4 -' for j = 0, 1, 2, • • • , 10, where p* max is 424, 1566, 

1027, 151, 100, 360, 162, 289, and 187 times larger than p, max at t = 0. 
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Fig. 20. The same as Fig. 4, but for non-adiabatic collapse for model (H) with K 1 ' 2 = 9.49. The 
contour lines are drawn for p*/p* max = lCP 0,4 ^ for j = 0, 1, 2, • • • , 10, where p* max is 165, 168, 
4276, 2684, 165, 73.6, 161, 1136, 1192, and 124 times larger than p, max at t = 0. The time 
appears in units of p^lL 2 at t = 0. 
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21. The same as Fig. 4, but for non-adiabatic collapse for model (K) with K~ x / 2 = 4.00. The 
contour lines are drawn for p*/p» max = 10~ ' 4 -' for j = 0, 1, 2, • • • , 10, where p, max is 45.9, 379, 
1366, 256, 23.1, and 219 times larger than p* max at t = 0. Note that p c = p max at t = in this 
case. 



